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We study the energy and the static spin structure factor of the ground state of the spin-1/2 
quantum Heisenberg antiferromagnetic model on the kagome lattice. By the iterative application of 
few Lanczos steps on accurate projected fermionic wave functions and the Green's function Monte 
Carlo technique, we find that a gapless (algebraic) U(l) Dirac spin liquid is competitive with previ- 
ously proposed gapped (topological) Z2 spin liquids. By performing a finite-size extrapolation of the 
ground-state energy, we obtain an energy per site E/J = —0.4365(2), which is equal, within three 
error-bars, to the estimates given by density-matrix renormalization group (DMRG). Our estimate 
is obtained for a translationally invariant system, and, therefore, does not suffer from boundary 
effects, like in DMRG. Moreover, on finite toric clusters at the pure variational level, our energies 
are lower compared to those from DMRG calculations. 

PACS numbers: 75.10.Jm, 75.10.Kt, 75.40.Mg, 75.50.Ee 



Introduction. The spin-1/2 quantum Heisenberg anti- 
ferromagnet (QHAF) on the kagome lattice provides a 
conducive environment to stabilize a quantum paramag- 
netic phase of matter down to zero temperature, [1-3] a 
fact that has been convincingly established theoretically 
from several studies, including exact diagonalization, [4— 
8] series expansion, [9, 10] quantum Monte Carlo, [11] 
and analytical techniques. [12] The question of the pre- 
cise nature of the spin-liquid state of the kagome spin- 
1/2 QHAF has been intensely debated on the theoreti- 
cal front, albeit without any definitive conclusions. Dif- 
ferent approximate numerical techniques have claimed 
a variety of ground states. On the one hand, density- 
matrix renormalization group (DMRG) calculations have 
claimed for a fully gapped (non-chiral) Z 2 topological 
spin-liquid ground state that does not break any point 
group symmetry. [13, 14] On the other hand, an algebraic 
and fully symmetric U(X) Dirac spin liquid has been pro- 
posed as the ground state, by using projected fermionic 
wave functions and variational Monte Carlo (VMC) ap- 
proach. [15-20] In addition, valence bond crystals have 
been suggested from many other techniques. In partic- 
ular, a 36-site unit cell valence-bond crystal [21-23] was 
proposed using quantum dimer models, [24-28] series ex- 
pansion [29, 30] and multi-scale entanglement renormal- 
ization Ansatz (MERA) [31] techniques. Finally, a recent 
coupled cluster method (CCM) suggested a q = (uni- 
form) state. [32] 

On general theoretical grounds, the Z 2 spin liq- 
uids in two spatial dimensions are known to be sta- 
ble phases, [47-49] as compared to algebraic 17(1) spin 
liquids, which are known to be only marginally sta- 
ble. [50] However, explicit numerical calculations using 
projected wave functions have shown the U(l) Dirac 
spin liquid to be stable (locally and globally) with re- 



spect to dimcrizing into all known valence-bond crys- 
tal phases. [15, 17, 18, 20] Furthermore, it was shown 
that, within this class of Gutzwiller projected wave func- 
tions, all the fully symmetric, gapped Z 2 spin liquids 
have a higher energy compared to the U{1) Dirac spin 
liquid. [19, 51] 

On the experimental front, the kagome spin-1/2 
QHAF model is well reproduced in Herbertsmithite 
(ZnCu3(OH) 6 Cl2), a compound with perfect kagome 
lattice geometry. [33-42] All experimental probes on 
Herbertsmithite point towards a spin-liquid behavior 
down to 20mK (i.e., four orders of magnitude smaller 
than the super-exchange coupling), which was estab- 
lished on the magnesium version of Herbertsmithite 
(MgCu3(OH) 6 Cl2). [43-45] Raman spectroscopic studies 
give further hints towards a gapless (algebraic) spin liq- 
uid. [46] 

In this letter, we systematically improve the projected 
fermionic wave functions of the U(l) Dirac and other 
competing spin liquids by applying a few Lanczos steps 
on large clusters, implemented stochastically within a 
variational Monte Carlo method. [52] We perform a zero- 
variance extrapolation of the energy and the static spin 
structure factor, which enables us to extract their exact 
values in the ground state on large cluster sizes and ob- 
tain an accurate estimate of the thermodynamic limit. In 
addition, we use Green's function Monte Carlo method, 
with the fixed-node (FN) approximation, [53] to extract 
the physical properties of the true ground state. Our 
main result is to show that the U(l) gapless spin liq- 
uid has an energy quite close to recent DMRG esti- 
mates, [13, 14] thus representing a very competitive state 
for the spin-1/2 QHAF on the kagome lattice (if not the 
true ground state). 

Model, wave functions, and numerical techniques. 
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FIG. 1. (Color online) Variational energies of the 48-site 
cluster as a function of the variance of the energy, for zero, 
one and two Lanczos steps. The ground-state energy is es- 
timated by extrapolating the three variational results to the 
zero- variance limit by a quadratic fit. Three different starting 
wave functions are used. The U(l) Dirac spin liquid has also 
been studied using FN approximation. 
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FIG. 2. (Color online) The same as in Fig. 1 for the 192-site 
cluster. Here, only Z2[0, 7r]/3 and U(l) Dirac states have been 
considered. 



The Hamiltonian for the spin-1/2 quantum Heisenberg 
antiferromagnetic model is: 



u 



(1) 



where J > and (ij) denotes sum over nearest-neighbor 
pairs of sites. The Si are spin-1/2 operators at each site 
i. All energies will be given in units of J. 

The physical variational wave functions are defined by 



FIG. 3. (Color online) The thermodynamic estimate of the 
ground-state energy obtained by a finite size extrapolation of 
the estimated ground-state energies (see Table I). The energy 
on the 36-site cluster is from exact diagonalization. Compar- 
ison is also made with recent DMRG estimates. [13, 14] 



projecting non-correlated fcrmionic states: 

|*VMcOftj,Ay,^,C)> =V G \^MF(Xij,^ij,IJ;0), ( 2 ) 

where Vg — rL(l~ n ',T n *^) * s ^ ne mu Gutzwiller projec- 
tor enforcing the one fermion per site constraint. Here, 
l^ / MF(Xiji Ay, /j,, £)} is the ground state of a mean-field 
Hamiltonian constructed out of Abrikosov fermions and 
containing hopping, chemical potential, and singlet pair- 
ing terms: 
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where a =t,-j-, Xij 
chemical potential [a, we also consider real and imagi- 
nary components of on-site pairing, which are absorbed 
in £. The spin-liquid phases are characterized by dif- 
ferent patterns of distribution of the underlying gauge 
fluxes through the plaquettes which are implemented by 
a certain distribution of the phases of Xij and Ajj on 
the lattice links, in addition one also needs to specify the 
on-site terms n and (. [48, 54] 

Here, we want to improve previous variational calcula- 
tions and approach, in a systematic way, the true ground 
state. This task can be achieved by the application of 
Lanczos steps: [52] 
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where the ctfc's are additional variational parameters. 
The convergence of |^ p _ls) to the exact ground state 
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l^cx) is guaranteed for large p provided the starting state 
is not orthogonal to |* cx ), i.e., for (^cxI^vmc) 7^ 0. 
However, on large cluster sizes, only a few steps can be 
efficiently performed and here we consider the case with 
p = 1 and p — 2 (jp = corresponds to the original 
starting variational wave function) . Subsequently, an es- 
timate of the exact ground-state energy may be achieved 
by the method of variance extrapolation: for sufficiently 
accurate states, we have that E w E cx + constant x a 2 , 
where E = (H)/N and a 2 = ((H 2 ) - (H) 2 )/N are the 
energy and variance per site, respectively. Whence, the 
exact ground-state energy E ex can be extracted by fitting 
E vs a 2 for p = 0, 1, and 2. 

The energy, its variance, and other physical proper- 
ties of the wave functions corresponding to p — 0, 1, 
and 2 Lanczos steps are obtained using standard VMC 
method. Moreover, the pure variational approach may 
be improved by using the FN approach, in which the 
high-energy components of the variational wave function 
are (partially) filtered out. [53] In particular, in the FN 
Monte Carlo method, the ground state of an auxiliary 
FN Hamiltonian is obtained and the approximation con- 
sists in assigning the nodal surface a priori, based upon a 
given guiding wave function, which is generally the best 
variational state. The energies obtained in this way are 
variational, [53] and hence we have a controlled approx- 
imation of the original problem. Here, the guiding wave 
function is obtained by optimizing the mean-field state 
of Eq. (2) using the method described in Refs. [55, 56]. 
Then, we find the best Lanczos parameters a p and finally 
we perform the VMC and FN Monte Carlo calculations 
for |*p_LS> with p = 0, 1, and 2. 

Results. We performed our variational calculations on 
toric clusters with mixed periodic-antiperiodic bound- 
ary conditions on the mean-field Hamiltonian of Eq. (3), 
which ensure non-degenerate wave functions at half fill- 
ing. We first consider the 48-site cluster (i.e., 4x4x3). 
As our starting (p = 0) variational wave functions, we 
take three different spin liquids, namely, (i) the U(l) 
Dirac spin liquid, which has a Fermi surface consisting 
of two points. [15, 16] The structure of the wave function 
is such that 10% of the configurations |x) (in which elec- 
trons reside on different sites of the lattice with given spin 
along the z direction) have zero weight (i.e., (x\^>) = 0). 
(ii) The uniform RVB spin liquid, which consists of a 
large circular spinon Fermi surface, [17] and has 35% of 
the configurations with zero weight, (iii) The Z2[0,7r]/3 
spin liquid, which is fully gapped [57] and has a negligi- 
ble (0.001%) number of configurations with zero weight. 
The zero- weight configurations are not visited by the ran- 
dom walk in the variational Monte Carlo method. The 
effect of two Lanczos steps on these wave functions is 
shown in Fig. 1 (see also Table I for the actual val- 
ues of the energies of the U(l) Dirac state). Our esti- 
mate of the ground-state energy on the 48-site cluster 
is thus E/J = —0.437845(4), which is comparable with 
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FIG. 4. (Color online) Intensity plot of the static spin struc- 
ture factor S(q) on the 192-site cluster. 

the DMRG estimate on a torus. [14, 58] Also the best 
pure variational energies are comparable within the two 
methods, see Table I. We want to stress the fact that the 
extrapolated energy is the same (within error bars) upon 
starting from all three wave functions. This is mainly due 
to the fact that, on relatively small clusters, few Lanczos 
steps are enough to filter out the high-energy components 
of the initial wave function and get a good estimation of 
the ground-state energy. 

On larger sizes, the extrapolations of U(l) and 
Z 2 [0,7r]/3 states deviate, the former one giving a slightly 
lower extrapolation, see Fig. 2 for the 192-site cluster. 
This fact suggests that the actual ground state is better 
described by a gapless algebraic U(l) Dirac state, rather 
than a gapped topological Z 2 spin liquid, as reported 
by DMRG calculations. In the following, for obtaining 
the ground-state energies on larger clusters we used only 
the U(l) Dirac wave function as the starting variational 
state. In Table I, we report our best results on different 
clusters (see the Supplementary Material for plots of the 
variance extrapolations on 108- (see Fig. 6) and 144-site 
(see Fig. 7) clusters. Also see Fig. 8). We would like to 
emphasize that our best variational energy on a 108-site 
cluster is significantly lower compared to the correspond- 
ing DMRG one, see Table I. 

By using the ground-state energy estimates on different 
cluster sizes, we performed a finite size extrapolation, see 
Fig. 3. Our final estimate for the energy of the infinite 
two-dimensional system is: 

E 2 J>/J = -0.4365(2). (5) 
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Size O-LS 1-LS 2-LS O-LS + FN 1-LS + FN 2-LS + FN Var. DMRG Est. Ground State 



48 


-0.4293510(4) 


-0.4352562(3) 


-0.436712(1) 


-0.432130(2) 


-0.435834(3) 


-0.436942(2) -0.4366 


-0.437845(4) 


108 


-0.4287665(4) 


-0.4341032(5) 


-0.435787(3) 


-0.431507(1) 


-0.434823(2) 


-0.436072(1) -0.4316 


-0.437178(9) 


144 


-0.4286959(5) 


-0.4337616(4) 


-0.435515(4) 


-0.4314455(8) 


-0.434544(2) 


-0.435839(9) 


-0.43698(2) 


192 


-0.4286749(4) 


-0.4334481(5) 


-0.435255(4) 


-0.431437(2) 


-0.434325(4) 


-0.435633(8) 


-0.43674(3) 



TABLE I. Energies of the U(l) Dirac spin liquid with p = 0, 1, and 2 Lanczos steps on different cluster sizes obtained by 
variational and FN Monte Carlo are given. In the penultimate column, we report the best variational DMRG energies. [59] 
The ground-state energy of the spin-1/2 QHAF estimated by us using zero- variance extrapolation of VMC energy values on 
different cluster sizes is marked in bold. 



This estimate is slightly higher (see Fig. 3) compared 
to DMRG extrapolations of Refs. [13, 14]. However, an 
energy estimate which is slightly lower by only few error- 
bars does not necessarily mean it is more accurate. We 
would like to stress that the same value for the extrapo- 
lated energy is obtained by using the FN approach (see 
Fig. 9 in the Supplementary Material). Moreover, it is 
worth mentioning that our energies are obtained with 
a state that has all the symmetries of the lattice, while 
DMRG states are non-uniform (due to boundary effects) . 

Let us now move to the calculation of the spin-spin 
correlations, which is defined by: 

ij R 

where N is the total number of sites, i,j — 1, 2, and 3 
label the three sites in the unit cell, R defines the Bravais 
lattice, and (R) is the real space spin-spin correlation 
function. 

The U(l) Dirac spin liquid is characterized by a power 
law (~ l/^ 4 ) decay of real space, long distance spin-spin 
correlations. [16] Here, we study the evolution of its static 
spin structure factor S(q) on the 192-site cluster under 
the action of one and two Lanczos steps and zero- variance 
extrapolation. Our estimate of the ground-state S(q) is 
obtained by a zero- variance extrapolation, see Fig. 11 in 
the Supplementary Material. The corresponding inten- 
sity plot of the extrapolated S(q) is shown in Fig. 4. One 
can clearly see that at large q, the spectral weight is 
concentrated on the corners of the hexagon, not very dif- 
ferently from what is found in a recent DMRG study. [14] 
However, what really matters is the behavior of S(q) for 
small q (namely at long distance) . Although the applica- 
tion of few Lanczos steps may not be sufficient to change 
the long-distance properties (because the Hamiltonian is 
a local operator), our calculations show that S(q) at small 
q remains practically unchanged under the action of one 
or two Lanczos steps and the subsequent zero-variance 
extrapolation (see Figs. 10 and 12 in the Supplementary 
Material). 

Summary. In summary, our systematic numerical 
study shows that competitive variational wave functions 
based upon Abrikosov fermions may be obtained. In- 
deed, our estimation for the energy of a gapless (alge- 



braic) U(l) Dirac spin liquid is very close to the recent 
DMRG results, [13, 14] which supported a fully gapped 
Z2 topological spin-liquid ground state. In this respect, 
our results lend support to the view that the exotic al- 
gebraic spin liquid can in fact occur as a true ground 
state of the spin-1/2 QHAF on the kagome lattice. Very 
recently, other approximate approaches proposed alter- 
native ground states with broken symmetries. [60, 61] 

We would like to mention that a further improvement 
of the variational wave function would require an intro- 
duction of local monopole fluctuations over the static 
mean-field state of Eq. (3). On small system sizes, such 
fluctuations were shown to lower the energy of the system 
within the Schwinger boson approach. [62] However, on 
large clusters, it is extremely difficult to construct work- 
able wave functions with (even static) topological defects. 
It is worth mentioning that the possibility of another en- 
ergetically competing state entering the game remains 
open, this is a chiral Z2 topological spin liquid [63] which 
has been proposed as the ground state within a Schwinger 
boson mean-field theory, [64] but whose projected wave 
function study remains to be done on large clusters such 
as 48 sites so as to enable a comparison with the U(l) 
Dirac spin liquid. Finally, the projected wave functions 
can also be constructed for chiral valence-bond crystal 
phases and it would be interesting to study their ener- 
getics, especially in light of the fact that they have been 
proposed as a competing ground state using generalized 
quantum dimer models. [28] 
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Supplementary Material 




FIG. 5. A lattice of corner sharing triangles giving rise to two elementary plaquettes, the triangle and the hexagon. The 
kagome lattice is the most frustrated among the 11 Archimedean tilings possible in 2D and has a coordination number of 4. 




Variance of energy (o ) 



FIG. 6. (Color online) The same as in Fig. 1 of the main paper, for the 108-site cluster. Here, only the U(l) Dirac state is 
considered and has been studied using both variational and fixed-node Monte Carlo methods. 
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FIG. 7. (Color online) The same as in Fig. 1 of the main paper, for the 144-site cluster. Here, only the U(l) Dirac state is 
considered and has been studied using both variational and fixed-node Monte Carlo methods. 
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FIG. 8. (Color online) The variance extrapolations on all cluster sizes for different starting variational wave functions are 
shown together for comparison. The inset shows the magnification near zero-variance. Note: The magnitude of the pairing in 
the Z2[0, n]/3 SL Ansatz is fixed to unity. 
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FIG. 9. (Color online) For all four cluster sizes, the fixed-node Monte Carlo energies have been extrapolated to the zero of 
the difference -Bvmc — £fn- The extrapolated values are within error bars of the zero- variance extrapolated values given in 
Table. I of the main paper, this shows that the VMC estimates are exact. 
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FIG. 10. (Color online) The Brillouin zone and the available crystal momenta for the 48 and the 192 site clusters are marked. 
Also, the 13 inequivalent momenta for the 192 site cluster are labelled with numbers. 
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FIG. 11. (Color online) Quadratic fit of zero-a extrapolation of the S(q) for the 13 inequivalent momenta on the 192-site 
cluster. The starting wave function is the 17(1) Dirac spin liquid on which one and two Lanczos steps have been performed. 
The numbering of the momenta is according to the above given Fig. 6. 
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FIG. 12. (Color online) On the 192-site cluster, the S(q) are plotted for the q points lying on the path marked with arrows in 
the Brillouin zone, namely for the (7(1) Dirac spin liquid and the corresponding wave functions obtained by applying one, and 
two Lanczos steps on it, and also the extrapolated ground state. The error bars of S(q) are of the order of 10 -4 . 



